1.1 用TTree 结构存储事件数据

目的:学习如何将事件以TTree的格式存储到ROOT文件

内容:

  1. 靶上生成中子,$\gamma$,入射到长条型塑料闪烁体探测器的不同位置和深度。
  2. 入射粒子形成的光信号向探测器两端传输,形成时间信号tu和td。
  3. 将入射粒子的原始信息,以及探测器信息逐事件存入TTree中.
    • 入射粒子原始信息:粒子种类、能量,粒子在探测器上的入射位置、入射深度等
    • 探测器信息: 两端时间信息tu、td,能量沉积信息qu,qd
  4. 打开root文件,进行基本的分析。

中子飞行时间谱测量原理

探测器参数

探测器的几何设置见下图:

  • 探测器长度:2L, 忽略探测器宽度。
  • 探测器厚度:$\Delta D$
  • 光在闪烁体的传播速度: $v_{sc}$
  • 光在闪烁体的衰减常数:$\lambda$
  • 粒子从起始位置到探测器中心的最小距离:D
  • 粒子入射位置:x
  • 飞行距离:$d=\sqrt{D^2+x^2}$

飞行时间:以靶为起点

中子飞行时间:$TOF_n(ns)=\frac{72 \cdot d(m)}{\sqrt{E_n(MeV)}}$

$\gamma$ 飞行时间:$TOF_{\gamma}(ns)=3.333\cdot d(m)$

探测器两端时间信号(不考虑其他offset):

  • $t_u=TOF+(L-x)/v_{sc}$
  • $t_d=TOF+(L+x)/v_{sc}$
  • $x=(t_d-t_u)*v_{sc}/2$
  • $TOF=(t_u+t_d)/2 -L/v_{sc}$

探测器两端能量(x 处能量沉积$Q_0$):

  • $Q_u=Q_0 exp[-(L-x)/\lambda]$
  • $Q_d=Q_0 exp[-(L+x)/\lambda]$
  • $x=\frac{\lambda}{2}ln\frac{q_u}{q_d}$
  • $Q_0=e^{2L/\lambda}*Q_u*Q_d$

代码

// 运行方法:
// 将下列代码保存到tree.cc
// 在ROOT环境中运行:
// .x tree.cc
// 生成 tree.root 文件
void tree(){
//常量声明
  const Double_t D=500.;//cm, distance between target and the scin.(Center)
  const Double_t L=100.;//cm, half length of the scin.
  const Double_t dD=5.;//cm, thickness of the scin.
  const Double_t TRes=1.;//ns, time resolution(FWHM) of the scintillator.
  const Double_t Lambda=380.;//cm, attenuation lenght of the scin.
  const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin. 
  const Double_t Vsc=7.5;//ns/cm, speed of light in the scin.
  const Double_t En0=100;//MeV, average neutron energy
  const Double_t EnRes=50.;//MeV, energy spread of neutron(FWHM)
  const Double_t Eg0=1;//MeV, gamma energy  
  const Double_t Rg=0.3;//ratio of gamma,ratio of neutron 1-Rg 
 //1. 定义新ROOT文件,声明新的Tree 
  TFile *opf=new TFile("tree.root","recreate");//新文件tree.root,指针 *opf
  TTree *opt=new TTree("tree","tree structure");//新tree,指针 *opt
 //2. 声明在tree结构中定义需要的变量分支
  Double_t x;//入射位置
  Double_t e;//能量
  int pid;    //粒子种类,n:pid=1,g:pid=0
  Double_t tof, ctof;//TOF:粒子实际飞行时间,cTOF:计算得到的TOF
  Double_t tu, td;
  Double_t qu, qd;

  Double_t tu_off=5.5;//time offset,//PMT的度越时间+电缆上的传输时间
  Double_t td_off=20.4;//time offset

  //3. 将变量地址添加到tree结构中
    //第一个参数为变量名称,第二个为上面定义的变量地址,第三个为变量的类型说明,D表示Double_t。
  opt->Branch("x", &x, "x/D");
  opt->Branch("e", &e, "e/D");
  opt->Branch("tof", &tof, "tof/D");
  opt->Branch("ctof",&ctof,"ctof/D");
  opt->Branch("pid", &pid, "pid/I");
  opt->Branch("tu", &tu, "tu/D");
  opt->Branch("td", &td, "td/D");
  opt->Branch("qu", &qu, "qu/D"); 
  opt->Branch("qd", &qd, "qd/D");  

// histogram,ROOT文件中除了TTree结构外,还可存储histogram,graph等
  TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
  TRandom3 *gr=new TRandom3(0);//声明随机数

  //4. 循环,计算变量的值,逐事件往tree结构添加变量值。
  for(int i=0;i<100000;i++){
    x=gr->Uniform(-L, L);//均匀入射在中子探测器表面.
    Double_t Dr=D+gr->Uniform(-0.5,0.5)*dD;//粒子在探测器厚度范围内均匀产生光信号
    Double_t d=TMath::Sqrt(Dr*Dr+x*x);//m, flight path
    if(gr->Uniform() < Rg) { //判断为gamma入射
       pid=0;
       e=Eg0;
       tof=3.333*(d*0.01);
    }
    else {  //neutron
        pid=1;
        e=gr->Gaus(En0, EnRes/2.35); // neutron
        tof=72./TMath::Sqrt(e)*(d*0.01);//ns
    }
    tu=tof+(L-x)/Vsc+gr->Gaus(0,TRes/2.35)+tu_off;
    td=tof+(L+x)/Vsc+gr->Gaus(0,TRes/2.35)+td_off;
    ctof=(tu+td)/2.;//simplified calculation.
    hctof->Fill(ctof);
    Double_t q0=e*gr->Uniform();//energy of recoil proton in plas. 0-En
    qu=q0*TMath::Exp(-(L-x)/Lambda);
    qu=gr->Gaus(qu,qu*QRes/2.35);
    qd=q0*TMath::Exp(-(L+x)/Lambda);
    qd=gr->Gaus(qd,qd*QRes/2.35);
    opt->Fill();//5.将计算好的变量值填到Tree中
  }
  // 6.将数据写入root文件中
  hctof->Write();
  opt->Write();
  opf->Close();
}

练习:

  • 理解上述1.-6.步骤的逻辑关系。自行练习如何将数据已TTree格式存到ROOT文件中。
  • 理解如何在模拟中加入探测器分辨率信息。

查看ROOT文件

  • juypter运行tree->Draw()时,需要有网络链接才能正常显示。
In [1]:
%jsroot on
In [2]:
TFile *ipf=new TFile("tree.root");//root -l tree.root
ipf->ls()//ROOT 环境下 > .ls
TFile**		tree.root	
 TFile*		tree.root	
  KEY: TH1D	hctof;1	neutron time of flight
  KEY: TTree	tree;1	tree structure
In [3]:
//观察tree的结构
tree->Print()
******************************************************************************
*Tree    :tree      : tree structure                                         *
*Entries :   100000 : Total =         6822821 bytes  File  Size =    5757028 *
*        :          : Tree compression factor =   1.18                       *
******************************************************************************
*Br    0 :x         : x/D                                                    *
*Entries :   100000 : Total  Size=     802585 bytes  File Size  =     614392 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.30     *
*............................................................................*
*Br    1 :e         : e/D                                                    *
*Entries :   100000 : Total  Size=     802585 bytes  File Size  =     565529 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.42     *
*............................................................................*
*Br    2 :tof       : tof/D                                                  *
*Entries :   100000 : Total  Size=     802645 bytes  File Size  =     742154 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.08     *
*............................................................................*
*Br    3 :ctof      : ctof/D                                                 *
*Entries :   100000 : Total  Size=     802675 bytes  File Size  =     741644 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.08     *
*............................................................................*
*Br    4 :pid       : pid/I                                                  *
*Entries :   100000 : Total  Size=     401467 bytes  File Size  =      42299 *
*Baskets :       13 : Basket Size=      32000 bytes  Compression=   9.48     *
*............................................................................*
*Br    5 :tu        : tu/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     755702 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.06     *
*............................................................................*
*Br    6 :td        : td/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     753799 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.06     *
*............................................................................*
*Br    7 :qu        : qu/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     769549 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.04     *
*............................................................................*
*Br    8 :qd        : qd/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     769490 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.04     *
*............................................................................*

观察指定事件数据

In [15]:
tree->Show(0);//显示指定的事件数据 event=1
======> EVENT:0
 x               = 99.3225
 e               = 92.424
 tof             = 38.1916
 ctof            = 63.5575
 pid             = 1
 tu              = 42.3357
 td              = 84.7793
 qu              = 86.1756
 qd              = 49.5036
In [5]:
tree->Show(100);
======> EVENT:100
 x               = -12.1667
 e               = 84.9707
 tof             = 39.2464
 ctof            = 65.5997
 pid             = 1
 tu              = 59.8935
 td              = 71.3059
 qu              = 48.0101
 qd              = 49.8971
In [6]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("td-tu>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [7]:
tree->Draw("tof>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
In [8]:
TH1D *hh=(TH1D*) ipf->Get("hctof"); //在代码中得到文件内数据指针的方法,在ROOT环境下课直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();

思考题1:观察上两个图有什么不同,分析其原因

In [13]:
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("ctof:td-tu>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw();

检查参数之间的关系是否与物理条件一致

  • gamma的到达探测器的时间(飞行时间)与探测器的x位置无关?
In [14]:
tree->Draw("ctof:td-tu>>(2000,-20,50,50,40,45)","pid==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯?
In [11]:
tree->Draw("td-tu:x","","colz");//观察两个参量的关联
c1->Draw();
In [12]:
tree->Draw("log(qu/qd):x","","colz");
c1->Draw();
In [3]:
!jupyter nbconvert 1.1_create_tree.ipynb --to html
[NbConvertApp] Converting notebook 1.1_create_tree.ipynb to html


[NbConvertApp] Writing 5013808 bytes to 1.1_create_tree.html


In [ ]: